cd('..')

clear all
addpath('.\functionsRA')
addpath('.\input')
%load('.\output\results_est_101_US_quart_mid.mat')
load('.\input\dists_US.mat','years','ny','last_month')
ending = '_US_monthquart';



load('.\input\results_est_101_US_monthquart.mat','Qs','pinew','R2')
Qs_q = Qs;
%pis_q = pis{1,1};
lth0 = size(Qs,1);

f=3;
lth = size(Qs_q,1);
% graph with comparison of the models
% p = plot((1:lth0)',100*Qs(:,1),(2:f:lth0-last_month)',100*Qs_q(:,1),(2:f:lth0-last_month)',100*mean(Qs((2:f:lth0-last_month),1)).*ones(lth,1),'--',(1:f:lth0-last_month),100*mean(Qs_q(:,1)).*ones(lth,1),'--','LineWidth',2);
% p(1).Color = '#A2142F'; p(2).Color = '#0072BD';
% p(3).Color = '#A2142F'; p(4).Color = '#0072BD';
% title('Result of Minimisations for US')
% legend('Model 101 monthly update of quarterly','Model 101 quarterly','Location','northwest')
% axis([1 lth0 0 12.5])
% set(gca,'Xtick',(1:12:lth0),'Ytick',(0:1:12),'XtickLabel',years,'XTickLabelRotation',45);
% ax = gca; ax.FontSize = 16; grid on
% ax.YAxis.TickLabelFormat = '%g%%';
% set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
% print('-bestfit',['.\figures matlab\100s_comparison' ending '.pdf'],'-dpdf')

% obs: in the graph above sometimes it seems that the quarterly is doing better, which is not possible, but that's just because those are months in which the quarterly model wasn't used (months 1,3,4,6,etc). 
% The dotted lines are the means of the fit in which both models are used, and monthly is slightly better. To check this, plot plot((1:44),[Qs((2:f:lth0-3),1) Qs_q]) and check they are basicaly the same, ans also that
% Qs((2:f:lth0-3),1) - Qs_q is alaways negative

mdl = 1;
name = ['model_' num2str(mdl+100)];

% graph with probabilities for the A matrix over time, begin careful about which parameters are assumed to be constant or not
% plot((1:lth0),100*pinew{mdl,1},'LineWidth',2);
% legend('\pi_{L}','\pi_{H}','\pi_{nL,t}','\pi_{nH,t}','\pi_{nn,t}','\pi_{mr}','Location','northwest','NumColumns',2)
% title(['Probabilities for Model (' num2str(mdl+100) ')'])
% axis([1 lth0 0 52])
% set(gca,'Xtick',(1:12:lth0),'Ytick',(0:5:50),'XtickLabel',years,'XTickLabelRotation',45);
% ax = gca; ax.FontSize = 16; grid on
% ax.YAxis.TickLabelFormat = '%g%%';
% set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
% print('-bestfit',['.\figures matlab\' name '' ending '\' name '' ending '_probs.pdf'],'-dpdf')

% graph with the prob of average inflation from t+6~t+10 being a disaster
pi_lim = [(-1:5)'; Inf]; pi_bar = [-2; (-0.5:4.5)'; 6];
prob_avg_dis = zeros(lth0,2); prob_avg_dis04 = prob_avg_dis; prob_avg_dis04_1to5 = prob_avg_dis; prob_avg_dis_1to5 = prob_avg_dis;
load('.\input\dists_US.mat','infl')
infl = infl(:);

for i = 1:lth0
    bin = (pi_lim>=infl(i)); bt = find(bin,1);
    x = pinew{1,1}(:,i);
    A = Amatrix(x,mdl);
    g = avg_inf_6to10(A,bt,pi_bar,pi_lim);
    prob_avg_dis(i,:) = [g(1) g(end)];
    prob_avg_dis04(i,:) = [g(1)+g(2) g(end-1)+g(end)];
%! new
    g = avg_inf_1to5(A,bt,pi_bar,pi_lim);
    prob_avg_dis_1to5(i,:) = [g(1) g(end)];
    prob_avg_dis04_1to5(i,:) = [g(1)+g(2) g(end-1)+g(end)];
end
% plot((1:lth0),100*prob_avg_dis,'LineWidth',2)
% title(['Prob of Avg Infl t+6~10 is a disaster, US, Model(' num2str(mdl+100) ')'])
% legend('Deflation','Hyperinflation','Location','northeast')
% axis([1 lth0 0 11])
% set(gca,'Ytick',(0:1:11),'XtickLabel',years,'XTickLabelRotation',45,'Xtick',(1:12:lth0));
% ax = gca; ax.FontSize = 16; grid on
% ax.YAxis.TickLabelFormat = '%g%%';
% set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
% print('-bestfit',['.\figures matlab\' name '' ending '\' name '' ending '_prob_avgdis.pdf'],'-dpdf')
% 
% % graph with the prob of average inflation from t+6~t+10 being a disaster with thresholds of 0 and 4 instead of -1 and 5
% plot((1:lth0),100*prob_avg_dis04,'LineWidth',2)
% title(['Prob of Avg Infl t+6~10 is a disaster, US, Model(' num2str(mdl+100) ')'])
% legend('Deflation','Hyperinflation','Location','northwest')
% axis([1 lth0 0 32.5])
% set(gca,'Ytick',(0:2:32),'XtickLabel',years,'XTickLabelRotation',45,'Xtick',(1:12:lth0));
% ax = gca; ax.FontSize = 16; grid on
% ax.YAxis.TickLabelFormat = '%g%%';
% set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
% print('-bestfit',['.\figures matlab\' name '' ending '\' name '' ending '_prob_avgdis04.pdf'],'-dpdf')
% 
% % graph with the prob of average inflation from t+1~t+5 being a disaster with thresholds of 0 and 4
% plot((1:lth0),100*prob_avg_dis04_1to5,'LineWidth',2)
% title(['Prob of Avg Infl t+1~t+5 is a disaster, US, Model(' num2str(mdl+100) ')'])
% legend('Deflation','Hyperinflation','Location','northwest')
% axis([1 lth0 0 42])
% set(gca,'Xtick',(1:12:lth0),'Ytick',(0:4:44),'XtickLabel',years,'XTickLabelRotation',45);
% ax = gca; ax.FontSize = 16; grid on
% ax.YAxis.TickLabelFormat = '%g%%';
% set(gcf,'Units','centimeters','PaperSize',[17 12],'position',[0 0 17 12])
% print('-bestfit',['.\figures matlab\' name '' ending '\' name '' ending '_prob_avgdis04_1to5.pdf'],'-dpdf')
% 
% % graph on R2 and RMSE
% p = plot((1:lth0),100*R2(:,1),'-',(1:lth0),R2(:,1)*0+100*mean(R2(:,1)),'r--','LineWidth',2);
% p(1).Color = '#A2142F';
% p(2).Color = '#A2142F';
% title(['Model Fit, US, Model(' num2str(100+mdl) ')'])
% axis([1 lth0 0 100])
% set(gca,'Xtick',(1:12:ny*12),'Ytick',(0:10:100),'XtickLabel',years,'XTickLabelRotation',45);
% ax = gca; ax.FontSize = 16; grid on
% ax.YAxis.TickLabelFormat = '%g%%';
% 
% yyaxis right
% p = plot((1:lth0),100*Qs,'-',(1:lth0),Qs*0+100*mean(Qs),'--','LineWidth',2);
% p(1).Color = '#0072BD';
% p(2).Color = '#0072BD';
% axis([1 lth0 0 20]); ax = gca; ax.FontSize = 16; box on; grid on
% ax.YAxis(2).Color = [0 0 0]; ax.YAxis(2).TickLabelFormat = '%g%%';
% set(gcf,'Units','centimeters','PaperSize',[17 12])
% set(gca,'Ytick',(0:2:20),'Units','centimeters','position',[1.8 3 11.3 7.3]); % [left bottom width height])
% legend('R2, lhs','mean R2','RMSE, rhs','mean RMSE','Location','southoutside','NumColumns',4,'Units','centimeters','position',[1.45 0.3 12 0.7])
% print('-bestfit',['.\figures matlab\' name '' ending '\' name '' ending '_fit.pdf'],'-dpdf')
% 

% Export probabilities to excel - this is the file that latter is imporated into STATA

filename = '.\input\prob_avginfl_monthly.xlsx';
filesheet = 'US_model_101_month';

label1 = {"5y5y", "5y5y", "next_5y", "next_5y","5y5y", "5y5y", "next_5y", "next_5y"};
label2 = {"year", "month","5y5y_def_0_4", "5y5y_hyp_0_4", "next_5y_def_0_4", "next_5y_hyp_0_4","5y5y_def_m1_5", "5y5y_hyp_m1_5", "next_5y_def_m1_5", "next_5y_hyp_m1_5"};
prob = [prob_avg_dis04 prob_avg_dis04_1to5 prob_avg_dis prob_avg_dis_1to5];

date_year = [];
for i = 1 : ny 
    date_year = [date_year ; years(i,1).*ones(12,1)];
end 
date_year = date_year(1:size(prob,1),1);

date_month = [];
for i = 1 : ny 
    date_month = [date_month ; (1:12)'.*ones(12,1)];
end 
date_month = date_month(1:size(prob,1),1);

date = [date_year date_month];
output = [date, prob];

writecell(label1,filename,"Sheet", filesheet,"Range","C1");
writecell(label2,filename,"Sheet", filesheet,"Range","A2");
writematrix(output,filename,"Sheet", filesheet,"Range","A3");

% Exporting parameters to excel

filename = '.\input\parameters_mdl101month.xlsx';
filesheet = 'US_mdl101_month';

par =  [date cell2mat(pinew(1,1))'];
label3 = {"year", "month", "pi_L","pi_H","pi_nL","pi_nH","pi_nn","pi_mr"};

writecell(label3,filename,"Sheet", filesheet,"Range","A1");
writematrix(par,filename,"Sheet", filesheet,"Range","A2");
